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Introduction 

Both  experimental  techniques  and  numerical  modeling  can  be  applied  to  study  the  same 
hydrodynamic  processes.  As  was  noted  in  [1],  the  development  and  verification  of  the 
diagnostics  methods  often  prove  to  be  the  most  laborious  part  of  the  numerical  experiment  as  in 
the  laboratory  experiment.  The  diagnostics  usually  involve  the  processing  of  the  results,  for 
example,  the  computation  of  integral  characteristics,  the  distribution  of  the  pressure  coefficient 
Cp  along  the  aircraft  wing  surfaces,  plotting  the  isolines,  etc.  These  diagnostics  are  made  with 
the  aid  of  corresponding  computer  codes,  which  enable  one  to  relate  the  discrete  variables  of 
modeling  to  the  quantities,  which  are  easier  to  interpret  or  compare  with  theory  or  experiment. 

The  problems  related  to  the  experimental  diagnostics  prove  to  be  more  difficult  than  in  the 
case  of  the  diagnostics  of  the  numerical  modeling  results.  The  numerical  experiment  diagnostics 
can  indeed  be  carried  out  without  any  connection  with  the  basic  computations,  whereas  the 
laboratory  measurements  usually  introduce  the  disturbances  in  experiment.  The  probes  can  alter 
the  local  flow,  absorb  the  heat  from  the  medium.  On  the  other  hand,  in  the  numerical 
experiment,  there  is  a  problem  of  the  tyranny  of  numbers:  what  to  do  with  the  dozens  or 
hundreds  of  thousands  of  numbers  obtained  at  the  output  of  the  numerical  experiment? 

The  localization  of  various  singularities  in  the  flow:  shock  waves,  slip  lines,  etc.  is  one  of 
the  problems  of  the  interpretation  of  numerical  results  in  aerodynamics.  Previously  a  theory  was 
presented  in  [2,  3],  which  enables  one  to  assess  the  accuracy  of  the  shock  wave  localization,  for 
example,  by  maximum  of  the  flow  velocity  gradient.  It  turned  out  that  the  divergence  property 
of  an  employed  finite-difference  scheme  is  one  of  the  necessary  conditions  for  ensuring  a  small 
error  of  various  gradient  techniques  for  shock  localization,  or  the  differential  analyzers,  as  they 
were  termed  in  [2,  3].  Then  it  is  possible  to  prove,  under  certain  additional  conditions,  the 
existence  and  uniqueness  of  such  a  point  in  the  zone  of  a  smeared  shock  wave,  the  finite- 
difference  shock  wave  center,  whose  location  coincides  with  the  exact  location  of  the  shock 
wave  front.  In  the  case  of  non-divergence  difference  schemes,  the  numerical  experiments  in  [1  - 
4]  have  shown  that  that  the  error  in  shock  wave  front  localization  may  be  significant  and  even 
be  a  quantity  of  the  order  0(1). 

Thus,  a  theoretical  base  was  formulated  in  [2,  3]  for  constructing  such  algorithms  for 
shock  localization  in  numerical  solutions,  which  possess  the  convergence,  that  is  their  error 
tends  to  zero  as  the  grid  steps  are  refined. 

Another  problem,  which  often  emerges  at  the  investigation  of  the  shockwave  structure  of 
flows  on  the  basis  of  numerical  experiment,  is  the  absence  of  a  priori  information  on  the 
location  and  types  of  singularities  in  the  spatial  computational  region.  To  solve  this  problem  it 
was  proposed  in  [5,  6]  for  the  first  time  to  apply  the  methods  of  digital  image  processing  for  the 
detection  of  strong  discontinuities  in  the  numerical  solutions  of  fluid  dynamics  problems.  In  [7, 
3]  these  methods  were  augmented  by  the  methods  of  digital  pattern  recognition  in  order  to 
classify  the  detected  singularities  lines  in  types.  The  performance  of  the  developed  algorithms 
for  the  localization  and  classification  of  singularities  was  demonstrated  in  [3]  on  a  number  of 
two-dimensional  aerodynamics  problems:  the  transonic  flow  around  the  airfoil,  the  double 
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Mach  reflection  of  a  shock  wave  on  a  wedge,  the  supersonic  flow  in  a  channel  with  a  lower  wall 
backward-facing  step.  The  pattern  recognition  methods  described  in  [3]  were  applied  in  [8]  for 
the  diagnostics  of  the  discontinuities  structure  in  the  problem  of  the  shock  wave  diffraction 
around  a  planar  wedge. 

It  should  be  noted  that  the  methods  of  digital  pattern  recognition  were  applied  in  [3,  5-8] 
only  for  the  diagnostics  of  pure  gas  flows.  At  the  same  time,  as  is  known  from  the  theory  of 
two-phase  flows  of  the  mixtures  of  gas  and  solid  particles,  the  shock  waves  as  well  as  other 
singularities  can  arise  in  such  flows,  and  the  number  of  different  types  of  singularities  in  two- 
phase  flows  even  proved  to  be  larger  than  in  the  case  of  pure  gas  flows  [9,  10]. 

A  generalization  of  the  algorithms  for  digital  localization  and  classification  of 
singularities  described  in  [2,  3,5-8]  is  proposed  in  the  present  work  for  the  case  of  two-phase 
flows  of  the  mixtures  of  gas  and  solid  particles.  The  examples  of  the  diagnostics  of  specific 
flows  are  presented. 

The  Problem  of  Shock  Wave  Interaction  with  a  Cloud  of  Particles 


Following  [11]  let  us  consider  a  rarefied  cloud  of  solid  spherical  particles  occupying  the 
region  Q2  (Fig.  1)  at  the  initial  moment  of  time  t  =  0.  At  some  time  t  >0,  a  shock  wave  (SW) 

impinges  from  the  left  on  the  cloud  of  particles.  In 
Q\  the  gas  flow  is  governed  by  the  Euler 
equations,  and  in  Q2  by  the  equations  of 
discrete/continual  model.  Choosing  the  x  axis 
along  the  SW  propagation,  and  the  y  axis  along  a 
normal  to  the  x  axis  let  us  write  a  system  of  the 
equations  for  a  flow  of  a  gas  -  particle  mixture  in 
x  planar  case: 
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Here  the  subscripts  1  and  2  denote  the  gas  and  particle  parameters,  respectively;  /  is  the 
distribution  function;  vj,  W\,  pu.  Pi,  Th  p,  y,  ex,  m\  are  the  gas  velocities  along  the  v  and  y,  the 
true  density,  mean  density,  temperature,  pressure,  the  adiabatic  exponent,  the  specific  internal 
energy,  and  the  volume  concentration  of  gas;  n  is  the  number  concentration  of  particles;  v2,  h'2> 
ax,  av,  p22,  r,  T2,  m2  are  the  velocities  and  accelerations  along  the  v  and  y,  the  true  density, 
radius,  temperature,  and  volume  concentration  of  particles;  /i,  X,  cv,  c  are  the  viscosity,  heat 
conduction,  heat  capacity,  and  sound  velocity  in  gas;  cs  is  the  specific  heat  of  the  particles 
material;  Re,  Nu,  Pr,  and  M12  are  the  Reynolds,  Nusselt,  Prandtl,  and  Mach  numbers. 

The  system  of  equations  (1)  was  solved  numerically  on  a  desktop  computer.  The  gas 
equations  were  solved  on  an  Eulerian  grid  with  the  third  order  of  accuracy,  and  the  kinetic 
equation  was  solved  in  the  Lagrange  variables.  The  equations  of  paths  and  heat  transfer  of 
particles  were  integrated  analytically  at  each  time  step.  This  has  enabled  us  to  increase  the 
accuracy  of  the  computation  of  particles  and  eliminated  a  limitation  for  the  time  step  r  <p\  1 
d2/p.  The  solid  wall  condition  was  set  on  yx  for  gas,  and  on  y2,  yi  ,  and  y4  the  symmetry 
condition  was  assumed.  The  mirror  reflection  condition  was  specified  for  particles  on  yh  and 
the  condition  for  particles  absorption  was  specified  on  y2,  Xi ,  and  y4  (see  Fig.  1). 
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The  Algorithm  for  Localization  and  Classification 
of  Singularities  in  Two-Phase  Flow 

As  was  shown  in  [9],  in  the  case  when  the  gas  density  is  much  smaller  than  the  density  of 
particles  material,  that  is  p\\/pn  «1,  m2  «1,  then  one  can  neglect  the  variation  of  temperature 
and  velocity  of  particles  across  the  shock  wave  front.  Then  the  relations  are  obtained  on  a  shock 
wave  in  two-phase  flow,  which  coincide  with  the  Rankine  —  Hugoniot  conditions  in  pure  gas: 

A+«„  -  D)  =  Pi  (V|'„  -£>); 

(2) 

Pi  (v\n  -  Df  +  Pi  <  =  A  (vln  -  D)2  +  A  »/,  ; 

A+  +Pi  1 P11  +  («i in  Jru\S~  !2  =  £\  +Pi  / Pu  +  (ui„  +uu)/2~ 

The  superscript  “+”  refers  here  to  the  medium  state  ahead  of  the  discontinuity  front,  and  the 
superscript  ”  refers  to  the  state  behind  the  discontinuity  front.  Since  the  relations  p\\/p22  1, 

«2  1  are  satisfied  in  the  case  of  the  problem  of  the  interaction  of  a  shock  wave  in  gas  with  the 

coal  dust  cloud,  one  can  use  for  the  localization  and  classification  of  strong  discontinuities  in 
the  given  two-phase  flow  the  algorithm,  which  was  developed  previously  in  [3,  7]  for  the  case 
of  pure  gas  flows.  In  this  connection,  we  describe  below  only  briefly  the  basic  stages  of  this 
algorithm.  These  are  the  following  two  stages: 

Stage  I.  Localization  of  lines  of  all  discontinuities. 

Stage  II.  The  classification  of  the  detected  discontinuity  lines  in  types  with  the  aid  of  the 
test  of  the  jump  relations  (2). 

The  procedure  for  localization  of  discontinuity  lines  consists  of  the  following  stages: 

1 .  Computation  of  the  gradient  magnitude  gy  in  all  nodes  (i,j)  of  the  digital  image  for  gas 
density  pUJ.  For  this  purpose  the  formulas  of  the  Sobel’s  edge  detector  [3,  5-7]  were  used: 

Sij  H  (Pij-ij+i  +2p\jj+ 1  +  A,/+i ,j+i ) -  ( Pi,i-\,j-\  +2p\jj-i  +  A, ;+i, y-i)  I 
+  1  (A  ,i+i,j+i  +  2-Pij+ij  +  Pij+ij-i )  (A,i-ij+i  +  2  Pij-ij  +  A,m,./-i  )  I  • 

1.  Thresholding. 

2.  Edge  thinning. 

3.  Tracing  contour  segments. 

As  regards  the  algorithm  for  contour  tracing,  it  was  not  yet  realized  in  [3].  This  algorithm 
was  implemented  in  [12]  to  detect  the  boundaries  of  the  multiply  connected  stability  regions  of 
difference  schemes.  Since  the  algorithm  for  tracing  contour  segments  described  in  [12]  is 
universal,  we  have  applied  it  without  changes  for  tracing  the  segments  of  discontinuity  lines  on 
the  basis  of  the  points  of  these  discontinuities,  which  were  detected  at  the  foregoing  stages.  The 
above  algorithm  is  based  on  the  ideas  for  tracing  contour  segments,  which  were  presented  in 
[13,14]. 

The  classification  of  the  lines  found  at  stage  1  into  types  of  discontinuities  was  performed 
at  stage  II  with  the  aid  of  the  sequential  classification  algorithm  described  in  [3].  In  this 
algorithm,  all  the  conditions  (2)  on  a  discontinuity  were  verified  sequentially.  To  identify  the 
shock  wave  fronts  also  the  condition  for  the  entropy  increase  behind  the  shock  wave  front  was 
checked.  The  algebraic  relations  valid  on  purely  contact  discontinuities  and  tangential 
discontinuities  were  checked,  which  are  satisfied  on  the  discontinuities  of  these  types,  cf.  [3]. 
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Fig.  2.  The  isolines  of  gas  density  p\. 


Fig.  3.  The  results  of  localization  and 
classification  of  singularities  in  the  gas 
phase  of  the  two-phase  flow:  solid  line  is 
the  normal  shock  wave,  dashed  line  is  the 
oblique  shock  wave,  the  dash-dot  line  is 
the  tangential  discontinuity,  the  black 
circles  show  a  purely  contact  discontinuity. 

1  2  3  4  5 

The  computation  of  two-phase  flow  and  subsequent  localization  and  classification  of 
singularities  were  executed  with  the  aid  of  the  corresponding  Fortran  programs.  Then  the 
plotting  of  isolines  of  gas  density  (see  Fig.  2)  and  the  lines  of  various  discontinuities  (see  Fig.  3) 
was  carried  out  with  the  aid  of  a  program  written  in  the  language  of  the  software  system 
Mathematica  4.1.  Figure  2  shows  the  isolines  of  the  gas  density  pi  at  the  moment  of  time  when 
the  shock  wave  has  already  passed  through  the  cloud  of  particles.  The  gas  parameters  ahead  of 
the  SW  correspond  to /jM0=l .3-10  3  g/cm3,  T®  =  293  K, p°  =  1  atm,  the  SW  Mach  number M0  = 
3.0.  The  initial  volume  concentration  of  particles  in  the  dust  cloud  at  t  =  0  is  m2  =  0.01.  The 
density  of  coal  particles  p22  =1.2  g/cm3,  the  particles  diameter  d  =  100  pm.  Figure  2  shows  40 
different  isolines  of  the  gas  density.  It  can  be  seen  that  the  isolines  coalesce  in  the  region  of  the 
impinging  shock  wave  and  the  shock  reflected  from  the  particles  cloud.  The  coalescence  of 
isolines  in  the  region  of  a  contact  discontinuity  behind  the  cloud  of  particles  is  not  so 
pronounced.  This  is  related  to  the  following  two  circumstances.  First,  a  relatively  crude  grid  of 
56x21  cells  was  used  in  the  computations  whose  results  are  presented  in  Figs.  2  and  3.  Second, 
the  contact  discontinuities  are  smeared  more  intensely  than  the  shock  waves  while  using  the 
shock-capturing  schemes  [2,  3]. 

Figure  3  presents  the  results  of  the  application  of  the  above  described  algorithm  for 
localization  and  classification  of  singularities  in  two-phase  flows  on  the  basis  of  shock¬ 
capturing  computational  results.  It  can  be  seen  that  the  normal  impinging  shock  wave  becomes 
slightly  curved  after  its  passage  through  the  cloud  of  particles.  A  contact  discontinuity  may  be 
seen  in  gas  between  the  cloud  of  particles  and  the  shock  wave.  It  is  a  purely  contact 
discontinuity  near  the  x  axis  and  it  goes  over  to  a  tangential  discontinuity  as  the  distance  from 
the  x  axis  increases.  This  tangential  discontinuity  separates  the  gas  percolated  through  the  cloud 
from  the  gas,  which  was  located  to  the  right  of  the  cloud  at  the  initial  moment  of  time  t  =  0  and 
then  accelerated  by  the  shock  wave. 

The  cloud  of  particles  was  plotted  in  Fig.  3  as  follows.  N0  marker  particles  ( N0>  1)  were 
assigned  to  the  initial  volume  concentration  of  particles  m2,  so  that  the  (1/  No)  th  fraction  of  the 
volume  concentration  of  particles  corresponds  to  each  marker  particle.  Then,  knowing  the 

volume  concentration  of  particles  m2j  in  each  cell  (i,  j),  we  can  find  the  number  of 
corresponding  marker  particles  tty  in  this  cell  by  formula  %  =  [m2jjN()  /  m2  ] ,  where  [a]  is  the 

integral  part  of  the  number  a ,  which  does  not  exceed  a.  The  coordinates  ( Xyhy^k )  of  marker 
particles  in  the  (i,  j)  cell  were  found  thereafter  by  formulas:  Xyy  =  (i  -  1  +  R)  hh  yjk  = 
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(j  -  1+R)h2,  where  R  is  a  pseudorandom  quantity  lying  in  the  interval  [0,  1].  This  quantity  was 
generated  in  the  Mathematicaa  4.1  program  with  the  aid  of  the  built-in  function  Random[  ].  It 
can  be  seen  from  Fig.  3  that  the  cloud  of  particles  at  the  time  under  consideration  has  not  yet 
been  deformed  significantly  in  comparison  with  its  initial  shape. 
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